library(tidyverse)

df <- haven::read_dta(here::here("PSRM_Dataset_Grapevine Effect.dta"))
hezbolla <- haven::read_dta(here::here("Corstange raw data_Grapevine Effect.dta"))
endline <- haven::read_dta(here::here("PSRM_RawData_Grapevine Effect.dta"))


# Bar plot of treatment and control groups --------------------------------

# Figure 1
df %>% 
  mutate(country_lab = case_when(country == 1 ~ "Chad",
                             country == 2 ~ "Niger",
                             country == 3 ~ "Burkina Faso")) %>% 
  group_by(d11_exp, country, d11_dummy, country_lab) %>% 
  summarise(n = n()) %>% 
  filter(!is.na(d11_exp)) %>% 

ggplot(.) +
  geom_bar(aes(x = forcats::fct_reorder(country_lab, country), y = n, fill = factor(d11_exp)), 
           stat = "identity", color = "white", size = 1) +
  facet_wrap(~d11_dummy, labeller = as_labeller(c(`0` = "Control group\n(n = 1,939)",
                          `1` = "Treatment group\n(n = 1,931)"))) +
  scale_fill_viridis_d(guide = guide_legend(title = ""), 
                       labels = c("Not at all", "Somewhat", "A lot", "Don't know", "Refused")) +
  theme_bw() +
  theme(panel.grid = element_blank()) +
  labs(y = "Observations", x = "")

ggsave(here::here("fig1_bar.png"),
       dpi = 800, height = 4.5, width = 9)


# Endorse 1 marginal effects ----------------------------------------------

# Figure 2

# Marginal effects estimated and exported from Stata. 

endorse1_margins <- 
  tibble::tribble(
         ~group,      ~est,       ~se,       ~low,     ~high, ~days,
      "Control", 0.2840428, 0.0309885,  0.2223968, 0.3456889,    0L,
    "Treatment",  0.579745, 0.0312693,  0.5175404, 0.6419495,    0L,
      "Control", 0.2655128, 0.0126399,   0.240368, 0.2906577,    5L,
    "Treatment", 0.5308956, 0.0131257,  0.5047844, 0.5570067,    5L,
      "Control", 0.2469828, 0.0286927,   0.189904, 0.3040617,   10L,
    "Treatment", 0.4820462, 0.0347983,  0.4128213,  0.551271,   10L,
      "Control", 0.2284528, 0.0543212,  0.1203906,  0.336515,   15L,
    "Treatment", 0.4331967, 0.0639818,  0.3059165,  0.560477,   15L,
      "Control", 0.2099228, 0.0808847,  0.0490173, 0.3708283,   20L,
    "Treatment", 0.3843473, 0.0939178,  0.1975149, 0.5711798,   20L,
      "Control", 0.1913928, 0.1076936, -0.0228441, 0.4056296,   25L,
    "Treatment", 0.3354979, 0.1240627,  0.0886977, 0.5822982,   25L,
      "Control", 0.1728627, 0.1346013,  -0.094902, 0.4406275,   30L,
    "Treatment", 0.2866485, 0.1542941, -0.0202916, 0.5935887,   30L,
      "Control", 0.1543327, 0.1615584, -0.1670584, 0.4757238,   35L,
    "Treatment", 0.2377991, 0.1845695, -0.1293684, 0.6049667,   35L,
      "Control", 0.1358027, 0.1885438, -0.2392709, 0.5108763,   40L,
    "Treatment", 0.1889497, 0.2148703, -0.2384958, 0.6163953,   40L
    )


ggplot(data = endorse1_margins) +
  geom_hline(yintercept = 0, color = "darkred", size = 1.5, alpha = 0.45) +
  geom_ribbon(aes(x = days, ymin = low, ymax = high, group = group), 
            size = 1.5, alpha = 0.2) +
  geom_line(aes(x = days, y = est, color = group), 
            size = 1.5) +
  geom_text(aes(x = 0, y = 0.5, label = "Treatment"),
            size = 4, hjust = 0, vjust = 1, check_overlap = T) +
  geom_text(aes(x = 0, y = 0.21, label = "Control"),
            size = 4, hjust = 0, vjust = 1, check_overlap = T) +
  geom_rug(data = filter(df, !is.na(endorse1_expbd) & !is.na(endorse1_dummyb)), 
                         aes(x = date_wave2, y = NA_real_),
           alpha = 0.25, position = position_jitter(width = 0.15)) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno",
                     guide = F, begin = 0, end = 0.5) +
  theme(panel.grid = element_blank()) +
  labs(x = "Days since first interview",
       y = "Prediction")

ggsave(here::here("fig2_endorse1.png"),
       dpi = 800, height = 5, width = 8)


# Married -----------------------------------------------------------------

# Figure 3
# Marginal effects estimated and exported from Stata. 

married1 <- 
  tibble::tribble(
    ~est,       ~se,       ~low,      ~high, ~days,
    -0.0724152, 0.0338654,  -0.141238, -0.0035923,    0L,
    -0.0032159, 0.0612406, -0.1276718,    0.12124,    5L,
    0.0659834, 0.1322292, -0.2027386,  0.3347055,   10L,
    0.1351827, 0.2057895, -0.2830318,  0.5533972,   15L,
    0.204382, 0.2799011, -0.3644455,  0.7732096,   20L,
    0.2735814, 0.3542183, -0.4462768,  0.9934395,   25L,
    0.3427807,  0.428634, -0.5283084,    1.21387,   30L
  )



m1 <- ggplot(data = married1) +
  geom_hline(yintercept = 0, color = "darkred", size = 1.5, alpha = 0.45) +
  geom_ribbon(aes(x = days, ymin = low, ymax = high), 
              size = 1.5, alpha = 0.2) +
  geom_line(aes(x = days, y = est), 
            size = 1.5) +
  geom_rug(data = filter(df, !is.na(endorse1_expd) & !is.na(ethmarry)), 
           aes(x = date_wave1, y = NA_real_),
           alpha = 0.25, position = position_jitter(width = 0.15)) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno",
                        guide = F, begin = 0, end = 0.5) +
  theme(panel.grid = element_blank()) +
  labs(x = "Days since first interview",
       y = "Marginal effect",
       title = "a) Survey 1")


married2 <-
  tibble::tribble(
          ~est,       ~se,       ~low,     ~high, ~days,
    -0.0418755,  0.042233, -0.1258903, 0.0421393,    0L,
    -0.0526889, 0.0362807, -0.1248627, 0.0194849,    5L,
    -0.0635022, 0.0429624,  -0.148968, 0.0219636,   10L,
    -0.0743155, 0.0580694, -0.1898341,  0.041203,   15L,
    -0.0851289,   0.07678, -0.2378687,  0.067611,   20L,
    -0.0959422, 0.0970314, -0.2889686, 0.0970841,   25L,
    -0.1067556, 0.1180332, -0.3415611,   0.12805,   30L,
    -0.1175689, 0.1394467, -0.3949728,  0.159835,   35L,
    -0.1283823, 0.1611079, -0.4488772, 0.1921126,   40L
    )

m2 <- ggplot(data = married2) +
  geom_hline(yintercept = 0, color = "darkred", size = 1.5, alpha = 0.45) +
  geom_ribbon(aes(x = days, ymin = low, ymax = high), 
              size = 1.5, alpha = 0.2) +
  geom_line(aes(x = days, y = est), 
            size = 1.5) +
  geom_rug(data = filter(df, !is.na(endorse1_expbd) & !is.na(ethmarry)), 
           aes(x = date_wave2, y = NA_real_),
           alpha = 0.25, position = position_jitter(width = 0.15)) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno",
                        guide = F, begin = 0, end = 0.5) +
  theme(panel.grid = element_blank()) +
  labs(x = "Days since first interview",
       y = "Marginal effect",
       title = "b) Survey 2")

ggsave(here::here("fig3_married.png"),
       dpi = 800, height = 5, width = 11,
       gridExtra::arrangeGrob(m1, m2, ncol = 2))




# Hezbolla ----------------------------------------------------------------

# Figure 4
# Marginal effects estimated and exported from Stata. 

hezbolla <- hezbolla %>% 
  mutate(hezbolla = case_when(q18 %in% c(-99, -98) ~ NA_real_, TRUE ~ q18,
                              q21 %in% c(-99, -98) ~ NA_real_, TRUE ~ q21,
                              q22 %in% c(-99, -98) ~ NA_real_, TRUE ~ q22,
                              q24 %in% c(-99, -98) ~ NA_real_, TRUE ~ q24)) %>% 
  mutate(date = case_when(q1 == "6/29/2011 0:00" ~ 0,
                          q1 == "6/30/2011 0:00" ~ 1,
                          q1 == "7/1/2011 0:00" ~ 2,
                          q1 == "7/2/2011 0:00" ~ 3,
                          q1 == "7/3/2011 0:00" ~ 4,
                          q1 == "7/4/2011 0:00" ~ 5,
                          q1 == "7/5/2011 0:00" ~ 6,
                          q1 == "7/6/2011 0:00" ~ 7,
                          q1 == "7/7/2011 0:00" ~ 8,
                          q1 == "7/8/2011 0:00" ~ 9,
                          q1 == "7/9/2011 0:00" ~ 10,
                          q1 == "7/10/2011 0:00" ~ 11,
                          q1 == "7/11/2011 0:00" ~ 12,
                          q1 == "7/14/2011 0:00" ~ 15),
         date2 = case_when(q2 == 1 ~ date - 2,
                           q2 == 2 ~ date - 0,
                           q2 == 3 ~ date - 0,
                           q2 == 4 ~ date - 1,
                           q2 == 5 ~ date - 1,
                           q2 == 6 ~ date - 0,
                           q2 == 7 ~ date - 12,
                           q2 == 8 ~ date - 0,
                           q2 == 9 ~ date - 3,
                           q2 == 10 ~ date - 6,
                           q2 == 13 ~ date - 8,
                           q2 == 15 ~ date - 7,
                           q2 == 16 ~ date - 7,
                           q2 == 17 ~ date - 9,
                           q2 == 18 ~ date - 9,
                           q2 == 19 ~ date - 9,
                           q2 == 20 ~ date - 11))

hezbolla <-
  hezbolla %>% 
  mutate(shianeigh = case_when(q30 == 3 ~ "Majority Shia", 
                               TRUE ~ "Minority Shia"))

# Figure 4
hez <- 
  tibble::tribble(
    ~group,     ~est,       ~se,      ~low,    ~high, ~days,
    "Minority Shia", 3.110172, 0.1541778,  2.777091, 3.443253,    1L,
    "Majority Shia", 3.307689, 0.2079774,  2.858381, 3.756997,    1L,
    "Minority Shia", 2.938333, 0.0671383,  2.793289, 3.083376,    4L,
    "Majority Shia", 3.728398, 0.2614311,   3.16361, 4.293185,    4L,
    "Minority Shia", 2.766493, 0.2613693,  2.201839, 3.331147,    7L,
    "Majority Shia", 4.149107, 0.3941608,  3.297574, 5.000639,    7L,
    "Minority Shia", 2.594653, 0.4638857,  1.592489, 3.596817,   10L,
    "Majority Shia", 4.569816, 0.5516377,  3.378075, 5.761556,   10L,
  )


ggplot(data = hez) +
  geom_hline(yintercept = 0, color = "darkred", size = 1.5, alpha = 0.45) +
  geom_ribbon(aes(x = days, ymin = low, ymax = high, group = group), 
              size = 1.5, alpha = 0.2) +
  geom_line(aes(x = days, y = est, color = group), 
            size = 1.5) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno", begin = 0, end = 0.5, guide = F) +
  geom_text(aes(x = 1, y = 2.75, label = "Minority Shia"),
            size = 4, hjust = 0, vjust = 1, check_overlap = T) +
  geom_text(aes(x = 1, y = 4.25, label = "Majority Shia"),
            size = 4, hjust = 0, vjust = 1, check_overlap = T) +
  geom_rug(data = filter(hezbolla, !is.na(q18)), 
           aes(x = date2, y = NA_real_),
           alpha = 0.25, position = position_jitter(width = 0.15)) +
  theme(panel.grid = element_blank()) +
  labs(x = "Days since first interview",
       y = "Prediction") +
  scale_x_continuous(limits  = c(1,10),
                     breaks = seq(1,10,2))

ggsave(here::here("fig4_hezbolla.png"),
       dpi = 800, height = 5, width = 7)



# Figure 5
hez2 <- 
  tibble::tribble(
    ~treatcontrol,          ~group,     ~est,       ~se,     ~low,    ~high, ~days,
        "control", "Minority Shia", 3.077938, 0.2031905, 2.638971, 3.516904,    1L,
        "control", "Majority Shia", 4.178338, 0.0490046,  4.07247, 4.284206,    1L,
        "control", "Minority Shia", 3.073309, 0.0111053, 3.049318, 3.097301,    4L,
        "control", "Majority Shia",  4.27347, 0.0310196, 4.206456, 4.340483,    4L,
        "control", "Minority Shia",  3.06868, 0.2209958, 2.591248, 3.546113,    7L,
        "control", "Majority Shia", 4.368601, 0.0457688, 4.269724, 4.467479,    7L,
        "control", "Minority Shia", 3.064052, 0.4329362,  2.12875, 3.999354,   10L,
        "control", "Majority Shia", 4.463733, 0.0761521, 4.299216,  4.62825,   10L,
      "treatment", "Minority Shia", 2.368303, 0.1447464, 2.055598, 2.681009,    1L,
      "treatment", "Majority Shia", 4.341362, 0.2147845, 3.877348, 4.805376,    1L,
      "treatment", "Minority Shia", 2.487531,  0.028942, 2.425006, 2.550056,    4L,
      "treatment", "Majority Shia", 4.527758, 0.2236347, 4.044625, 5.010891,    4L,
      "treatment", "Minority Shia", 2.606759, 0.1451965, 2.293081, 2.920437,    7L,
      "treatment", "Majority Shia", 4.714154, 0.2395907,  4.19655, 5.231758,    7L,
      "treatment", "Minority Shia", 2.725986, 0.2858054, 2.108541, 3.343431,   10L,
      "treatment", "Majority Shia",  4.90055, 0.2613544, 4.335928, 5.465172,   10L,
    )


ggplot(data = hez2) +
  geom_hline(yintercept = 0, color = "darkred", size = 1.5, alpha = 0.45) +
  geom_ribbon(aes(x = days, ymin = low, ymax = high, group = treatcontrol), 
              size = 1.5, alpha = 0.2) +
  geom_line(aes(x = days, y = est, color = treatcontrol), 
            size = 1.5) +
  geom_rug(data = filter(hezbolla, !is.na(q18)), 
           aes(x = date2, y = NA_real_, group = shianeigh),
           alpha = 0.25, position = position_jitter(width = 0.15)) +
  theme_bw() +
  scale_color_viridis_d(option = "inferno", begin = 0, end = 0.5, guide = guide_legend(title = "")) +
  theme(panel.grid = element_blank(),
        legend.position = "top") +
  labs(x = "Days since first interview",
       y = "Prediction") +
  scale_x_continuous(limits  = c(1,10),
                     breaks = seq(1,10,2)) +
  facet_wrap(~group)


ggsave(here::here("fig5_hezbolla.png"),
       dpi = 800, height = 5, width = 7)





# Endline raw bar ---------------------------------------------------------

endline %>% 
  mutate(country_lab = case_when(id2 == 1 ~ "Chad",
                                 id2 == 2 ~ "Niger",
                                 id2 == 3 ~ "Burkina Faso"),
         
         d11_exp = case_when(is.na(b_d11traitement) ~ b_d11controle,
                             TRUE ~ b_d11traitement),
         d11_dummy = case_when(!is.na(b_d11controle) ~ 0,
                               !is.na(b_d11traitement) ~ 1)) %>% 
  group_by(id2, d11_exp, country_lab, d11_dummy, country_lab) %>% 
  summarise(n = n()) %>% 
  filter(!is.na(d11_dummy)) %>%
  
  ggplot(.) +
  geom_bar(aes(x = forcats::fct_reorder(country_lab, id2), y = n, fill = factor(d11_exp)), 
           stat = "identity", color = "white", size = 1) +
  facet_wrap(~d11_dummy, labeller = as_labeller(c(`0` = "Control group\n(n = 1,963)",
                                                  `1` = "Treatment group\n(n = 1,954)"))) +
  scale_fill_viridis_d(guide = guide_legend(title = ""), 
                       labels = c("Not at all", "Somewhat", "A lot", "Don't know", "Refused")) +
  theme_bw() +
  theme(panel.grid = element_blank()) +
  labs(y = "Observations", x = "") +
  scale_y_continuous(limits = c(0, 800))

ggsave(here::here("figA2_bar.png"),
       dpi = 800, height = 4.5, width = 9)
